Spatial Interpolation
Matthew Whittle
16/09/2021
Packages
library(raster)
library(tidyverse)
library(sf)
library(gstat)
library(sp)
library(lubridate)
library(getarc)
library(janitor)
library(vroom)
library(fs)
library(gridExtra)
library(leaflet)
library(leafem)
library(conflicted)
conflict_prefer("select", "dplyr")
conflict_prefer("filter", "dplyr")Loading the data
The data is split into 6 csv files, one for each region in the dataset. I’m loading them by listing the files, using the vroom function to read and combine them, then converting the names to ‘snake_case’.
Note: the vroom function loads a representation of the CSVs and doesn’t load the whole csv into memory until it is called by R. This is much faster than using
read_csv.
temps <-
# List files
fs::dir_ls("project_data/", regexp = "daily*") %>%
# Load them into a dataframe
vroom::vroom(show_col_types = FALSE) %>%
# tidy the column names
janitor::clean_names()
# Preview
head(temps)Data Quality
The data has some quality issues to address. There are some suspect values in the data such as temperatures below -50 and years before 1995. I’m going to fix them with an assumption or convert them to missing values.
set.seed(1)
year_hist <-
temps %>%
slice_sample(n = 10000) %>%
ggplot(aes(x = year)) +
geom_histogram() +
labs(title = "Year") +
theme_minimal()
temp_hist <-
temps %>%
slice_sample(n = 10000) %>%
ggplot(aes(x = avg_temperature)) +
geom_histogram() +
labs(title = "Avg Temp") +
theme_minimal()
grid.arrange(year_hist, temp_hist)
I’m going to assume that days with a value of 0 should be 1; any years before 1995 are errors and should be missing values; and temperatures below -50 should also be missing values. Here I’m using the if_else function to replace suspect values with missing values.
I could solve this problem using interpolation too but I think that is probably beyond the scope of the analysis.
temps <-
temps %>%
mutate(day = if_else(day < 1, 1, day),
year = if_else(year < 1995, NA_real_, year),
avg_temperature = if_else(avg_temperature < -50, NA_real_, avg_temperature)
)There are some duplicate values on some days so I am summarising the temperatures here to remove them. This is important because I want to count the days with dangerous temperatures.
temps <-
temps %>%
mutate(date = ymd(glue::glue("{year} {month} {day}"))) %>%
group_by(region, country, city, year, month, day, date) %>%
slice(1)
head(temps)City Locations
I want to model the temperature changes across regions, to do this I’m going to add in a spatial component to the dataset. I’ve ge0code the locations of each city using another which I’ll load, turn into a simple feature and bind with the data in a future step.
city_xy <- read_csv("project_data/location/city-xy.csv")
city_points <- city_xy %>% st_as_sf(coords = c("longitude", "latitude"), crs = 4326)
city_points <- city_points %>% st_transform(4326)
world_boundaries <- st_read("project_data/location/world-boundaries.geojson")## Reading layer `world-boundaries' from data source
## `/Users/matthew/Documents/Data Analysis/NCW2021-MW/project_data/location/world-boundaries.geojson'
## using driver `GeoJSON'
## Simple feature collection with 249 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -89 xmax: 180 ymax: 83.6236
## Geodetic CRS: WGS 84leaflet() %>%
addProviderTiles(providers$OpenStreetMap) %>%
addFeatures(world_boundaries, popup = ~country) %>%
addFeatures(city_points, popup = ~city, color = "pink")Humidity Data
The risk level posed by different temperatures is dependent on the relative humidity (measured: 0-100%). I’ve used some regional monthly average humidity averages to develop some rasters. I’m loading the data as a stars object then getting the nearest value to each city.
humidity_rasters<- brick("project_data/humidity-monthly-average.grd")
plot(humidity_rasters$March)
I’m doing it this way instead of extracting the values because many of the cities are coastal and therefore return NA for the humidity values (as the sea is NA in my rasters).
humidity_points <-
stars::read_stars("project_data/humidity-monthly-average.grd") %>% st_as_sf()
city_humidity <-
bind_cols(city_points,
humidity_points[st_nearest_feature(city_points, humidity_points),] %>% st_drop_geometry())
city_humidity <-
city_humidity %>% st_drop_geometry() %>%
pivot_longer(
cols = -c(city, region, country, address),
names_to = "month",
values_to = "average_humidity"
)
temp_hum <-
temps %>%
ungroup() %>%
mutate(month_label = month(date, label = TRUE, abbr = FALSE)) %>%
left_join(city_humidity,
by = c("city", "region", "country", "month_label" = "month"))
head(temp_hum)Heat Index
The heat index is an approximate index of ‘feels like’ temperatures. This is dependent on relative humidity and temperature. The index can be used to define approximate risk levels from feels like temperatures. More detail is available on this website.
temp_hum <-
temp_hum %>%
mutate(
heat_index =
weathermetrics::heat.index(
t = avg_temperature,
rh = average_humidity * 100,
temperature.metric = "celcius"
),
danger_level =
case_when(heat_index > 51.67 ~ "Extreme Danger",
heat_index > 39.44 ~ "Danger",
heat_index > 32.22 ~ "Extreme Caution",
heat_index > 26.67 ~ "Caution",
TRUE ~ "Low Risk") %>%
factor(levels = c(
"Low Risk",
"Caution",
"Extreme Caution",
"Danger",
"Extreme Danger"
), ordered = TRUE)
)
head(temp_hum)Changes in Dangerous Heat Days over time?
Assuming humidity levels have stayed consistent across the globe since 1995 (which is a big assumption), how many days per year has each city experienced in each heat category?
annual_risk_days <-
temp_hum %>%
group_by(year, region, country, city, danger_level) %>%
count(name = "days")
annual_risk_days <- city_points %>% right_join(annual_risk_days)
head(annual_risk_days)Kriging
risk_2019 <-
annual_risk_days %>%
filter(year == 2019 & danger_level == "Extreme Danger")
risk_1995 <-
annual_risk_days %>%
filter(year == 1995 & danger_level == "Extreme Danger")Making a grid
ep_continents <- "https://services.arcgis.com/P3ePLMYs2RVChkJx/arcgis/rest/services/World_Continents/FeatureServer/0"
continents <- query_layer(endpoint = ep_continents, in_geometry = st_union(annual_risk_days))
continents_region <- continents %>% st_filter(annual_risk_days)
world_grid <- st_read("project_data/location/world-grid.geojson")## Reading layer `world-grid' from data source
## `/Users/matthew/Documents/Data Analysis/NCW2021-MW/project_data/location/world-grid.geojson'
## using driver `GeoJSON'
## Simple feature collection with 1337 features and 0 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -178 ymin: -87 xmax: 178 ymax: 81
## Geodetic CRS: WGS 84world_grid <- world_grid %>% st_transform(crs = 4326)
world_grid <- world_grid %>% st_filter(continents_region)
plot(continents_region$geometry) 
Kriging
require(automap)
krige_to_raster <-
function(krige, crs = 4326) {
raster <-
krige$krige_output %>%
st_as_sf() %>%
st_transform(crs) %>%
add_coords() %>%
select(X = x, Y = y, Z = var1.pred) %>%
st_drop_geometry() %>%
as.data.frame() %>%
raster::rasterFromXYZ()
raster::crs(raster) <- sp::CRS(glue::glue("EPSG:{crs}"))
return(raster)
}
add_coords <-
function(x, names = c("x", "y")) {
xy <- sf::st_coordinates(x)
xy <- tibble::as_tibble(xy)
colnames(xy) <- names
dplyr::bind_cols(x, xy)
}Annual Risk Changes